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Abstract 

The mathematical theory of pattern formation in electrically coupled networks of excitable neurons 
forced by small noise is presented in this work. Using the Freidlin-Wentzell large deviation theory for 
randomly perturbed dynamical systems and the elements of the algebraic graph theory, we identify and 
analyze the main regimes in the network dynamics in terms of the key control parameters: excitability, 
coupling strength, and network topology. The analysis reveals the geometry of spontaneous dynamics 
in electrically coupled network. Specifically, we show that the location of the minima of a certain 
continuous function on the surface of the unit n— cube encodes the most likely activity patterns generated 
by the network. By studying how the minima of this function evolve under the variation of the coupling 
strength, we describe the principal transformations in the network dynamics. The minimization problem 
is also used for the quantitative description of the main dynamical regimes and transitions between 
them. In particular, for the weak and strong coupling regimes, we present asymptotic formulae for the 
network activity rate as a function of the coupling strength and the degree of the network. The variational 
analysis is complemented by the stability analysis of the synchronous state in the strong coupling regime. 
The stability estimates reveal the contribution of the network connectivity and the properties of the 
cycle subspace associated with the graph of the network to its synchronization properties. This work is 
motivated by the experimental and modeling studies of the ensemble of neurons in the Locus Coeruleus, 
a nucleus in the brainstem involved in the regulation of cognitive performance and behavior. 



1 Introduction 

Direct electrical coupling through gap-junctions is a common way of communication between neurons, as 
well as between cells of the heart, pancreas, and other physiological systems ifTOl . Electrical synapses are 
important for synchronization of the network activity, wave propagation, and pattern formation in neuronal 
networks. A prominent example of a gap-junctionally coupled network, whose dynamics is thought to be 
important for cognitive processing, is a group of neurons in the Locus Coeruleus (LC), a nucleus in the 
brainstem lHHHI^. Electrophysiological studies of the animals performing a visual discrimination test 
show that the rate and the pattern of activity of the LC network correlate with the cognitive performance 
[381 . Specifically, the periods of the high spontaneous activity correspond to the periods of poor perfor- 
mance, whereas the periods of low synchronized activity coincide with good performance. Based on the 
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physiological properties of the LC network, it was proposed that the transitions between the periods of 
high and low network activity are due to the variations in the strength of coupling between the LC neurons 
1.38 J . This hypothesis motivates the following dynamical problem: to study how the dynamics of electrically 
coupled networks depends on the coupling strength. This question is the focus of the present work. 

The dynamics of an electrically coupled network depends on the properties of the attractors of the lo- 
cal dynamical systems and the interactions between them. Following [38J, we assume that the individual 
neurons in the LC network are spontaneously active. Specifically, we model them with excitable dynamical 
systems forced by small noise. We show that depending on the strength of electrical coupling, there are 
three main regimes of the network dynamics: uncorrelated spontaneous firing (weak coupling), formation 
of clusters and waves (intermediate coupling), and synchrony (strong coupling). The qualitative features 
of these regimes are independent from the details of the models of the individual neurons and network 
topology. Using the center manifold reduction |[8l[T8l and the Freidlin-Wentzell large deviation theory |[T4l . 
we derive a variational problem, which provides a useful geometric interpretation for various patterns of 
spontaneous activity. Specifically, we show that the location of the minima of a certain continuous function 
on the surface of the unit n— cube encodes the most likely activity patterns generated by the network. By 
studying the evolution of the minima of this function under the variation of the control parameter (coupling 
strength), we identify the principal transformations in the network dynamics. The minimization problem is 
also used for the quantitative description of the main dynamical regimes and transitions between them. In 
particular, for the weak and strong coupling regimes, we present asymptotic formulae for the activity rate as 
a function of the coupling strength and the degree of the network. The variational analysis is complemented 
by the stability analysis of the synchronous state in the strong coupling regime. In analyzing various aspects 
of the network dynamics, we pay special attention to the role of the structural properties of the network in 
shaping its dynamics. We show that in weakly coupled networks, only very rough structural properties of 
the underlying graph matter, whereas in the strong coupling regime, the finer features, such as the algebraic 
connectivity and the properties of the cycle subspace associated with the graph of the network, become 
important. Therefore, this paper presents a comprehensive analysis of electrically coupled networks of ex- 
citable cells in the presence of noise. It complements the existing studies of related deterministic networks 
of electrically coupled oscillators (see, e.g., ||TTl|T9l|26l|33l and references therein). 

The outline of the paper is as follows. In Section [2} we formulate the biophysical model of the LC net- 
work. Section|3]presents numerical experiments elucidating the principal features of the network dynamics. 
In Section |4j we reformulate the problem in terms of the bifurcation properties of the local dynamical sys- 
tems and the properties of the linear coupling operator. We then introduce the variational problem, whose 
analysis explains the main dynamical regimes of the coupled system. In Section [5| we analyze the stability 
of the synchronous dynamics in the strong coupling regime, using fast-slow decomposition. The results of 
this work are summarized in Section [6l 



2 The model 

2.1 The single cell model 

According to the dynamical mechanism underlying action potential generation, conductance-based models 
of neurons are divided into Type I and Type II classes |35|. The former assumes that the model is near the 
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Figure 1: a) The phase plane for ( |2.1| ), ( |2.2| ): nullclines plotted for the deterministic model (a = 0) 
and a trajectory of the randomly perturbed system (a > 0). The trajectory spends most time in a small 
neighborhood of the stable fixed point. Occasionally, it leaves the basin of attraction of the fixed point to 
generate a spike, b) The voltage timeseries, v{t), corresponding to spontaneous dynamics shown in plot a. 



saddle-node bifurcation, while the latter is based on the Andronov-Hopf bifurcation. Electrophysiological 
recordings of the LC neurons exhibit features that are consistent with the Type I excitability. The existing 
biophysical models of LC neurons use Type I action potential generating mechanism In accord with 

these experimental and modeling studies, we use a generic Type I conductance-based model to simulate the 
dynamics of the individual LC neuron 

Cv = -Iion{v,n) ^ aw, (2.1) 
noo{v)-n 

n . (2.2) 

t{v) 

Here, dynamical variables v{t) and n{t) are the membrane potential and the activation of the potassium 
current, Ik, respectively. C stands for the membrane capacitance. The ionic currents lioni^^ ^) are modeled 
using the Hodgkin-Huxley formalism (see Appendix for the definitions of the functions and parameter values 



used in ( |2.1| ) and ( |2.2| )). A small Gaussian white noise is added to the right hand side of ( |2.1| ) to simulate 
random synaptic input and other possible fluctuations affecting system's dynamics. Without noise (a = 0), 
the system is in the excitable regime. For a > 0, it exhibits spontaneous spiking. The frequency of the 
spontaneous firing depends on the proximity of the deterministic system to the saddle-node bifurcation and 
on the noise intensity. A typical trajectory of ( |2.1| ) and ( |2.2[ ) stays in a small neighborhood of the stable 



equilibrium most of the time (Fig. [T^). Occasionally, it leaves the vicinity of the fixed point to make a 
large excursion in the phase plane and then returns to the neighborhood of the steady state (Fig.[T^). These 
dynamics generate a train of random spikes in the voltage time series (Fig.[T]3). 

In neuroscience, the (average) firing rate provides a convenient measure of activity of neural cells and 
neuronal populations. It is important to know how the firing rate depends on the parameters of the model. 
In this paper, we study the factors determining the rate of firing in electrically coupled network of neurons. 
However, before setting out to study the network dynamics, it is instructive to discuss the behavior of the 



single neuron model first. To this end, we use the center-manifold reduction to approximate ( |2.1| ) and ( |2.2[ ) 
by a li? system: 

i = -U\z) + awu U{z) = /iz - ^z^ + ^fi^/^ (2.3) 
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Figure 2: The numerical approximation of the density of the time between the successive spikes in the 
voltage timeseries v{t), obtained by integration (2.1) and (2.2). The interspike intervals are distributed 
approximately exponentially. 



where z{t) is the rescaled projection of {v{t)^n{t)) onto a ID slow manifold, > is the distance to the 
saddle-node bifurcation, and a > is the noise intensity after rescaling. We postpone the details of the 
center-manifold reduction until we analyze a more general network model in §4.1| 



The time between two successive spikes in voltage time series corresponds to the first time the trajectory 



of (2.3) with initial condition z(0) = < overcomes potential barrier U{^/J1) — U{—^/Jl). The large 
deviation estimates (cf. O) yield the logarithmic asymptotics of the first crossing time r 

lim InE = 2U{^) = ^ E x exp J I , (2.4) 

where E zq stands for the expected value with respect to the probability generated by the random process 
z{t) with initial condition z{0) = zq. Throughout this paper, we use x to denote logarithmic asymptotics. 
It is also known that the first exit time r is distributed exponentially as shown in Fig.[2](cf. l[T2ll ). 

Equation ( |2.4| ) implies that the statistics of spontaneous spiking of a single cell is determined by the 



distance of the neuronal model ( |2.1[ ) and ( |2.2| ) to the saddle-node bifurcation and the intensity of noise. 
Below we show that, in addition to these two parameters, the strength and topology of coupling are important 
factors determining the firing rate of the coupled population. 



2.2 The electrically coupled network 



The network model includes n cells, whose intrinsic dynamics is defined by ( |2.1[ ) and ( |2.2| ), coupled by 
gap-junctions. The gap-junctional current that Cell i receives from the other cells in the network is given by 



= 1 



(t;y)-t;«), 



(2.5) 
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Figure 3: Examples of coupling schemes: (a) nearest neighbor (cf. Example [XT] ); (b) 2— nearest neighbor 
(cf. Examplel^l); (c) all-to-all (cf. Example [23]). 



where ^ > is the gap-junction conductance and 

1, Cell i and Cell j are connected, 



an = 0, (i, i) G In] . 
0, otherwise. ^ \ ^JJ l j 

Adjacency matrix A = (a^j) G M^^^ defines the network connectivity. By adding the coupling current to 
the right hand side of the voltage equation ( |2.1| ) and combining the equations for all neurons in the network, 
we arrive at the following model 

n 

= -Iio„(i;«,n«)+5^ai,(t;y)-t;«)+a«;«, (2.6) 



i=l 



n^(t;(^))-n(^) 
r(^;(*)) 



n« = , (2.7) 



where w^'^'^ are n independent copies of the standard Brownian motion. 



The network topology is an important parameter of the model ( |2.6| ) and ( |2.7| ). The following terminology 
and constructions from the algebraic graph theory |6| will be useful for studying the role of the network 
structure in shaping its dynamics. Let G = {V{G), ^{G)) denote the graph of interactions between the 
cells in the network. Here, V{G) = {vi, ^2, • • • , ^n} and E{G) = {ei, 62, . . . , e^} denote the sets of 
vertices (i.e., cells) and edges (i.e., the pairs of connected cells), respectively. Throughout this paper, we 
assume that G is a connected graph. For each edge ej = {^311^32) ^ ^{G) x V{G), we declare one of 
the vertices Vj-^ , Vj^ to be the positive end (head) of e^, and the other to be the negative end (tail). Thus, we 
assign an orientation to each edge from its tail to its head. The coboundary matrix of G is defined as follows 
(cf. M) 

1, Vj is a positive end of e^, 
— 1, Vj is a negative end of e^, (2.8) 
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0, otherwise. 



Let G = (y{G)^ E{G)) C G be a spanning tree of G, i.e., a connected subgraph of G such that \V{G) \ = n, 
and there are no cycles in G |6|. Without loss of generality, we assume that 

£;(G) = {ei,e2,...,e,_i}. (2.9) 

Denote the coboundary matrix of G by 
Matrix 

L = H^H (2.10) 
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Figure 4: Regular versus random connectivity. Both graphs in (a) and (b) have degree 4. The graph in (a) 
is formed using regular coupling scheme, whereas edges of the graph in (b) are generated using a random 
algorithm (cf. Example |2.4[). 




is called a graph Laplacian of G. The Laplacian is independent of the choice of orientation of edges that 
was used in the definition of i7 Q. Alternatively, the Laplacian can be defined as 



where D = diag{deg(i;i), deg(i;2), . . . deg{vn)} is the degree map and A is the adjacency matrix of G. 



denote the eigenvalues of L arranged in the increasing order counting the multiplicity. The spectrum of the 
graph Laplacian captures many structural properties of the network (cf. O [71 [9l)- In particular, the first 
eigenvalue of L, Ai(L) = 0, is simple if and only if the graph is connected |[T6ll . The second eigenvalue 
a = MiL) is called the algebraic connectivity of G, because it yields a lower bound for the edge and 
the vertex connectivity of G |[T6ll . The algebraic connectivity is important for a variety of combinatorial, 
probabilistic, and dynamical aspects of the network analysis. In particular, it is used in the studies of the 
graph expansion p2l, random walks 0, and synchronization of dynamical networks |[T7l[29l . 

Next, we introduce several examples of the network connectivity including nearest neighbor arrays of 
varying degree and a pair of degree 4 symmetric and random graphs. These examples will be used to 
illustrate the role of the network topology in pattern formation. 

Example 2.1. The nearest-neighbor coupling scheme is an example of the local connectivity (Fig. [jj^J. For 
simplicity, we consider a ID array. For higher dimensional lattices, the nearest neighbor coupling is defined 
similarly. In this configuration, each cell in the interior of the array is coupled to two nearest neighbors. 
This leads to the following expression for the coupling current: 




L^D-A, 



(2.11) 



Let 



Ai(L)<A2(L)<...<A,(L) 



/O') = ^(z;0'+l) - + ^(i;(^'-l) - v^^^), j = 2, 3, . . . 



n — 1. 



The coupling currents for the cells on the boundary are given by 



/(I) =^(^;(2) _^(i)) and 
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The corresponding graph Laplacian is 



( 1 
-1 



-1 
2 -1 








\ 




(2.12) 



V ... -1 1 y 

Example 2.2. The k— nearest neighbor coupling scheme is a natural generalization of the previous example. 
Suppose each cell is coupled to k of its nearest neighbors from each side whenever they exist or as many as 
possible otherwise: 

mm{fc,n— j'} ram{k,j} 



J = 2,3, . . . ,n - 1, 



(2.13) 



i=l i=l 

where we use a customary convention that Yl'j=a ' ' ' 
from ( |2.i3| ). 

Example 2.3. The all-to-all coupling features global connectivity (Fig. ^): 



ifb < a. The coupling matrix can be easily derived 



1,2,3, ... ,n. 



(2.14) 



i=l 



The Laplacian in this case has the following form 

/ n - 1 -1 

-1 n - 1 



V 



-1 



-1 



n 



-1 \ 
-1 

-1 y 



(2.15) 



The graphs in the previous examples have different degrees: ranging from 2 in Example [2^1] to n — 1 in 
Example |2.3| In addition to the degree of the graph, the pattern of connectivity itself is important for the 
network dynamics. This motivates our next example. 

Example 2.4. Consider a pair of degree 4 graphs shown schematically in Fig. |?] The graph in Fig. ^ has 
symmetric connections. The edges of the graph in Fig. ^ were selected randomly. Both graphs have the 
same number of nodes and equal degrees. 

Graphs with random connections like the one in the last example represent expanders, a class of graphs 
used in many important applications in mathematics, computer science and other branches of science and 
technology (cf. [22J). In Section [5] we show that dynamical networks on expanders have very good synchro- 
nization properties (see also L29J ). 

Example 2.5. Let {Gn} be a family of graphs on n vertices, with the following property: 

HGn) >a>0, n G N. (2.16) 

Such graphs are called (spectral) expanders 122\ \36\l . There are known explicit constructions of expanders, 
including the celebrated Ramanujan graphs ^28l \27\l . In addition, families of random graphs have good 
expansion properties. In particular, it is known that 

Prob {a2(G^) >d- 2Vd-l - e} = 1 - o^(l) Ve > 0, (2.17) 

where Gn stands for the family of random graphs of degree d > 3 and 1 / [731/ . 
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Figure 5: The fundamental relation of the rate of spontaneous activity and the coupling strength. The 
graphs in (a) are plotted for three coupling configurations: the nearest neighbor (dashed line), the 2— nearest 
neighbor coupling (solid line) and all-to-all coupling (dash-dotted line) (see Examples |2.ip.3| ). The graphs 
in (b) are plotted for the symmetric and random degree 4 graphs in solid and dashed lines respectively (see 
Example |2.4| ). (c) The firing rate plot for the model, in which the coupling is turned off for values of the 
membrane potential above the firing threshold. The symmetric (solid line) and random (dashed line) degree 
4 graphs are used for the two plots in (c). 



3 Numerical experiments 



The four parameters controlling the dynamics of the biophysical model p.6| ) and ( |2.7| ) are the excitability, 
the noise intensity, the coupling strength, and the network topology. Assuming that the system is at a fixed 
distance from the bifurcation, we study the dynamics of the coupled system for sufficiently small noise 
intensity a. Therefore, the two remaining parameters are the coupling strength and the network topology. 
We focus on the impact of the coupling strength on the spontaneous dynamics first. At the end of this section, 
we discuss the role of the network topology. The numerical experiments of this section show that activity 
patterns generated by the network are effectively controlled by the variations of the coupling strength. 



3.1 Three phases of spontaneous activity 

To measure the activity of the network for different values of the control parameters, we will use the average 
firing rate - the number of spikes generated by the network per one neuron and per unit time. Fig. [5^ shows 
that the activity rate varies significantly with the coupling strength. The three intervals of monotonicity 
of the activity rate plot reflect three main stages in the network dynamics en route to complete synchrony: 
weakly correlated spontaneous spiking, formation of clusters and wave propagation, and synchronization. 
We discuss these regimes in more detail below. 

Weakly correlated spontaneous spiking. For ^ > sufficiently small, the activity retains the features 
of spontaneous spiking in the uncoupled population. Fig. [6]3 shows no significant correlations between the 
activity of distinct cells in the weakly coupled network. The distributions of the interspike intervals are 
exponential in both cases (see Fig.|6](c,d)). There is an important change, however: the rate of firing goes 
down for increasing values of ^ > for small g. This is clearly seen from the graphs in Fig. [5] The 
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Figure 6: Spontaneous activity in uncoupled (a) and weakly coupled (b) networks. The corresponding 
distributions for the time intervals between successive spikes are exponential (c,d) with a slightly heavier 
tail in the latter case. 
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Figure 7: Coherent structures in the weakly coupled network: a) clusters and short waves, b,c) robust waves, 
d) nearly synchronous discharge. Networks shown in Figures 6][8 are coupled through the nearest neighbor 
scheme. 
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Figure 8: Clusters generated by the modified model, in which the electrical current from a given cell is 
turned off once the cell has crossed the threshold (see text for details). These experiments show that Factor 
A vs. B is responsible for forming clusters in the weak coupling regime. 
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decreasing firing rate for very weak coupling can also be noted from the interspike interval distributions in 
Fig.[6]:,d: the density in Fig.[6jl has a heavier tail. Thus, weak electrical coupling has a pronounced inhibitory 
(shunting) effect on the network dynamics: it drains the current from a neuron developing a depolarizing 
potential and redistributes it among the cells connected to it. This effect is stronger for networks with 
greater number of connections. The three plots shown in Fig. |5^ correspond to nearest-neighbor coupling, 
2— nearest neighbor coupling, and all-to-all coupling. Note that the slope at zero is steeper for networks 
with greater degree. 

Coherent structures. For increasing values of ^ > the system develops clusters, short waves, and ro- 
bust waves (see Fig. [7]). The appearance of these spatio-temporal patterns starts in the middle of the first 
decreasing portion of the firing rate plot in Fig. |5^ and continues through the next (increasing) interval of 
monotonicity. While patterns in Fig.|7]feature progressively increasing role of coherence in the system's dy- 
namics, the dynamical mechanisms underlying cluster formation and wave propagation are distinct. Factors 
A and B below identify two dynamical principles underlying pattern formation in this regime. 

Factor A: At the moment when one neuron fires due to large deviations from the rest state, neurons con- 
nected to it are more likely to be closer to the threshold and, therefore, are more likely to fire within a short 
interval of time. 

Factor B: When a neuron fires, it supplies neurons connected to it with depolarizing current. If the coupling 
is sufficiently strong, the gap-junctional current triggers action potentials in these cells and the activity 
propagates through the network. 

Factor A follows from the variational interpretation of the spontaneous dynamics in weakly coupled 
networks, which we develop in Section |4j It is responsible for the formation of clusters and short waves, 
like those shown in Fig. |7^. To show numerically that that Factor A (vs. Factor B) is responsible for the 



formation of clusters, we modified the model (|2.6|) and (2.7) in the following way. Once a neuron in the 



network has crossed the threshold, we turn off the current that it sends to the other neurons in the network 



until it gets back close to the the stable fixed point. We will refer to this model as the modified model (2.6) 



and ( |2.7| ). Numerical results for the modified model in Fig.[8^,b, show that clusters are formed as the result 



of the subthreshold dynamics, i.e., are due to Factor A. Factor B becomes dominant for stronger coupling. 
It results in robust waves with constant speed of propagation. The mechanism of the wave propagation 
is essentially deterministic and is well known from the studies of waves in excitable systems (cf. 1241). 
However, in the presence of noise, the excitation and termination of waves become random (see Fig.[7|b,c)). 

Synchrony. The third interval of monotonicity in the graph of the firing rate vs. the coupling strength is 
decreasing (see Fig. [5^). It features synchronization, the final dynamical state of the network. In this regime, 
once one cell crosses the firing threshold the entire network fires in unison. The distinctive feature of this 
regime is a rapid decrease of the firing rate for increasing g (see Fig. [5^). The slowdown of firing in the 
strong coupling regime was studied in |(32l| (see also EIIEHIMI)- When the coupling is strong the effect 
of noise on the network dynamics is diminished by the dissipativity of the coupling operator. The reduced 
effect of noise results in the decrease of the firing rate. In §5.4[ we present analytical estimates characterizing 
denoising by electrical coupling for the present model. 
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Figure 9: Spontaneous activity patterns generated by regularly (a) and randomly (b) connected degree-4 
networks (cf. Example [2!4| ) for the same value of the coupling strength g = 0.006. The randomly connected 
network is already synchronized (b), while the regular network is en route to synchrony (a). 



3.2 The role of the network topology 

All connected networks of excitable elements (regardless of the connectivity pattern) undergo the three dy- 
namical regimes, which we identified above for weak, intermediate, and strong coupling. The topology 
becomes important for quantitative description of the activity patterns. In particular, the topology affects the 
boundaries between different phases. We first discuss the role of topology for the onset of synchronization. 
The transition to synchrony corresponds to the beginning of the third phase and can be approximately iden- 
tified with the location of the point of maximum on the firing rate plot (see Fig. |5^,b). The comparison of 
the plots for 1— and 2— nearest-neighbor coupling schemes shows that the onset of synchrony takes place 
at a smaller value of g for the latter network. This illustrates a general trend: networks with greater number 
of connections tend to have better synchronization properties. However, the degree is not the only structural 
property of the graph that affects synchronization. The connectivity pattern is important as well. Fig. [9] 
shows that a randomly connected degree 4 network synchronizes faster than its symmetric counterpart (cf. 



Example 2.4). The analysis in ^4.4 shows that the point of transition to synchrony can be estimated using 
the algebraic connectivity of the graph a. Specifically, the network is synchronized, if 7 > a~^, where 7 
stands for the coupling strength in the rescaled nondimensional model. The algebraic connectivity is easy 



to compute numerically. For many graphs with symmetries including those in Examples |2.ip.3[ the al- 



gebraic connectivity is known analytically. On the other hand, there are effective asymptotic estimates of 
the algebraic connectivity available for certain classes of graphs that are important in applications, such as 
random graphs lITSl and expanders |[22ll . The algebraic connectivities of the graphs in Examples |2.1||Z2 



a = 0(n~^) tend to zero as n ^ oc. Therefore, for such networks one needs to increase the strength of 
coupling significantly to maintain synchrony in networks growing in size. This situation is typical for sym- 
metric or almost symmetric graphs. In contrast, it is known that for the random graph from Example |24] the 
algebraic connectivity is bounded away from zero (with high probability) as n ^ 00 ||T5l|22l. Therefore, 
one can guarantee synchronization in dynamical networks on such graphs using finite coupling strength 
when the size of the network grows without bound. This counter-intuitive property is intrinsic to networks 
on expanders, sparse well connected graphs |[22l[36l. For a more detailed discussion of the role of network 
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topology in synchronization, we refer the interested reader to Section 5 in ll29ll . 



The discussion in the previous paragraph suggests that connectivity is important in the strong coupling 
regime. It is interesting that to a large extent the dynamics in the weak coupling regime remains unaffected 
by the connectivity. For instance, the firing rate plots for the random and symmetric degree-4 networks 
(Example [5]) shown in Fig. |5]3 coincide over an interval in g near 0. Furthermore, the plots for the same 



pair of networks based on the modified model ( |2.6| ) and ( |2.7| ) are almost identical, regardless the disparate 
connectivity patterns underlying these networks. The variational analysis in §4.3| shows that, in the weak 
coupling regime, to leading order the firing rate of the network depends only on the number of connections 
between cells. The role of the connectivity in shaping network dynamics increases in the strong coupling 
regime. 



4 The variational analysis of spontaneous dynamics 



In this section, we analyze dynamical regimes of the coupled system ( |2.6| ) and ( |2.7| ) under the variation of 
the coupling strength. In §4.1[ we derive an approximate model using the center manifold reduction. In 
§4.2[ we relate the activity patterns of the coupled system to the minima of a certain continuous function 
on the surface of an n— cube. The analysis of the minimization problem for weak, strong, and intermediate 
coupling is used to characterize the dynamics of the coupled system in these regimes. 



4.1 The center manifold reduction 



In preparation for the analysis of the coupled system (2.6) and (2.7), we approximate it by a simpler system 



using the center manifold reduction |[8l [TH. To this end, we first review the bifurcation structure of the 
model. Denote the equations governing the deterministic dynamics of a single neuron by 

x = f(x,/x), (4.1) 

where x G and f : x ^ is a smooth function and is a small parameter, which controls the 
distance of (14. ID from the saddle-node bifurcation. 



Assumption 4.1. Suppose that at ji — {), the unperturbed problem \4.1\ has a nonhyperbolic equilibrium at 
the origin such that i?f (0, 0) has a single zero eigenvalue and the rest of the spectrum lies to the left of the 
imaginary axis. Suppose further that at ji — {) there is a homoclinic orbit to O entering the origin along the 
ID center manifold. 

Then under appropriate nondegeneracy and transversality conditions on the local saddle-node bifurcation 
at /X = 0, for ji near zero the homoclinic orbit is transformed into either a unique asymptotically stable 



periodic orbit or to a closed invariant curve having two equilibria: a node and a saddle [8| (Fig. 10). 
Without loss of generality, we assume that the latter case is realized for small positive fi, and the periodic 
orbit exists for negative fi. Let /i > be a sufficiently small fixed number, i.e., (|4.1|) is in the excitable 



regime (Fig. 10). For simplicity, we assume that the stable node near the origin is the only attractor of (4.1 ). 
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Figure 10: Local system (4.1 ) is near the saddle-node on an invariant circle bifurcation. 



We are now in a position to formulate our assumptions on the coupled system. Consider n local systems 



(4.1 ) that are placed at the nodes of the connected graph G = {V{G)^ E{G)), \V{G)\ = n, and coupled 
electrically: 

X = F(X, ^) - g{L J)X + a{In P)W, (4.2) 

where X = (x^^) , x^^) , . . . , x(^) )^ G x • • • x = R^^, F(X, /x) = (f (x^^) , /x) , f (x^^) , ^) , . . . , f (x ^ , /x)) ^, 

is an n X n identity matrix, P G R^^^, and L G R^^^ is the Laplacian of G. Matrix J G R^^^ defines 
the linear combination of the local variables engaged in coupling. In the the neuronal network model above, 
J = diag(l, 0). Parameters g and a control the coupling strength and the noise intensity respectively. W 
is a Gaussian white noise process in R^^. The local systems are taken to be identical for simplicity. The 
analysis can be extended to cover nonhomogeneous networks. 



We next turn to the center manifold reduction of ( |4.2| ). Consider ( |4.2| )o (the zero subscript refers to 
(J = 0) for /i = ^ = 0. By our assumptions on the local system ( |4.1[ ), i?f (0, 0) has a ID kernel. Denote 

e G ker Df{0, 0) /{O} and p G ker(i:>f (0, 0))"^ such that p^e^l. (4.3) 



By the center manifold theorem, there is a neighborhood of the origin in the phase space of ( |4.2[ ), B, and 
6 > such that for < S and \g\ < S, in B, there exists an attracting locally invariant n— dimensional 
slow manifold M^^g. The trajectories that remain in B for sufficiently long time can be approximated by 



those lying in M/^^g. Thus, the dynamics of (4.2 )o can be reduced to M/^^g, whose dimension is d times 



smaller than that of the phase space of (4.2 )o. The center manifold reduction is standard. Its justification 



relies on the Lyapunov-Schmidt method and Taylor expansions (cf. lUl). Formally, the reduced system is 



obtained by projecting ( |4.2[ )o onto the center subspace of ( |4.2[ )o for = ^ = (see [25]): 

y = aiy'^ - a2fi - a^gLy + 0(|yp, ^^), (4.4) 

where y = (yi, ^2, • • • , ^n) ^ R^, := (y^, • • • ? 1/^)5 provided that the following nondegeneracy 
conditions hold 

1 



ai = ^Q^P^H^e,0)\u=o ^0. (4.5) 
-|^j9Tf(0,/x) 1^.0 7^0, (4.6) 



a2 

as = p^Je^O. (4.7) 
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Conditions (|4.5|) and (|4.6|) are the nondegeneracy and transversality conditions of the saddle-node bifurcation 



in the local system ( |4.1| ). Condition ( |4.7[ ) guarantees that the projection of the coupling onto the center 
subspace is not trivial. All conditions are open. Without loss of generality, assume that nonzero coefficients 
<^i,2,3 are positive. 

Next, we include the random perturbation in the reduced model. Note that near the saddle-node bifur- 



cation (0 < /i <C 1), the vector field of (4.2 )o is much stronger in the directions transverse to M/j^^g than 
in the tangential directions. The results of the geometric theory of randomly perturbed fast-slow systems 



imply that the trajectories of (4.2) with small positive a that start close to the node of (4.2 )o remain in a 
small neighborhood of A^^,^ on finite intervals of time with overwhelming probability (see |3 | for specific 
estimates). To obtain the leading order approximation of the stochastic system ( |4.2| ) near the slow manifold, 
we project the random perturbation onto the center subspace of (4.2)o for = ^ = and add the resultant 
term to the reduced equation (|4.4[): 



y = aiy' - a2fi - a^gLy + aBW + . . . , B = 1^ ® {p' P) e M^^^^. (4.8) 

We replace BW by identically distributed a^w, where if; is a white noise process in and = \P'^p\' 
Here, | • | stands for the Euclidean norm of P^p G M^. After rescaling the resultant equation and ignoring 
the higher order terms, we arrive at the following reduced model 



(4.9) 



where w stands for a standard Brownian motion in and In = (1,1,... ,1) G M^. Here, with a slight 
abuse of notation, we continue to use a to denote the small parameter in the rescaled system. In the remain- 



der of this paper, we analyze the reduced model (4.9 ). 



4.2 The exit problem 



In this subsection, the problem of identifying most likely dynamical patterns generated by ( |4.2[ ) is reduced 
to a minimization problem for a smooth function on the surface of the unit cube. 



Consider the initial value problem for ( |4.9| ) 

i = f (z) - -fLz + aw, L = H^H, z{0) = zq e D cR"", (4.10) 

where 

i{z) = {f{zi)J{z2),...J{Zn)), f{0=e-h (4.11) 

and 

D = {z={zi,Z2,...,Zn):-2-b<Zi<l,ie [n] := {1,2,..., n}}, (4.12) 
where auxiliary parameter b > will be specified later. Let 

d+D = {z = (zi, Z2, . . . , Zn) : zeDk{3ie [n] = 1)} , (4.13) 

denote a subset of the boundary of D, dD. If z{t) G d^D, then at least one of the neurons in the network is 
at the firing threshold. It will be shown below that the trajectories of (4.10) exit from D through d^D with 
probability 1 as a ^ 0, provided 6 > is sufficiently large|^ Therefore, the statistics of the first exit time 

T = inf {t > : z{t) G dD} (4.14) 



Positive parameter b in the definition of D (cf. |4.12l) is used to exclude the possibihty of exit from D through dD \d^D. 
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and the distribution of the location of the points of exit z{t) E d^D characterize the statistics of the 
interspike intervals and the most probable firing patterns of ( |2.1| ) and ( |2.2| ), respectively. The Freidlin- 
Wentzell theory of large deviations [14J yields the asymptotics of r and z{t) for small a > 0. 



z = - — U^{z) + awt, (4.15) 



To apply the large deviation estimates to the problem at hand, we rewrite ( |4.10[ ) as a randomly perturbed 
gradient system 

d_ 
dz^ 
where 

U^{z) = ^{Hz, Hz) + $(z), $(z) = ^(^^)' ^ 3 + ^ " 3^'' ^^'^^^ 

i=l 

where (•, •) stands for the inner product in R^~^. The additive constant 2/3 in the definition of the potential 
function F(^) is used to normalize the value of the potential at the local minimum F(— 1) = 0. 



The following theorem summarizes the implications of the large deviation theory for ( |4.15 ) 



Theorem 4.2. Let z^^\ z^^^, . . . , z^^\ /c G N, denote the points of minima ofU^{z) on dD 

U^{z^^) = U min U^{z), i = 1, 2, . . . , A:, 

zedD 

and Z — UiLiI^^^^}- Then for any zq ^ D and 6 > 0, 

A) lim^^o^ Mz(r),Z) < 5} = 1, (4.17) 

B) lim^^o ^2 InE = U, (4.18) 

C) lim^^o F zo {exp{a-2(f7 - h)} < r < exp{a-^{U + h)}} = 1, V/i > 0, (4.19) 

where p(-, •) stands for the distance in W^. 

The statements A)-C) can be shown by adopting the proofs of Theorems 2.1, 3.1, and 4.1 of Chapter 4 
of O to the case of the action functional with multiple minima. 



Theorem |42] reduces the exit problem for ( |4.10| ) to the minimization problem 

U^{z) min, z G dD. (4.20) 
In the remainder of this section, we study ( |4.20[ ) for the weak, strong, and intermediate coupling strength. 

4.3 The weak coupling regime 

In this subsection, we study the minima of U^{z) on dD for small |7|. First, we locate the points of minima 



small I7I (cf. Theorem 4.4). 



of the U^{z) for 7 = (cf. Lemma 4.3 ). Then, using the Implicit Function Theorem, we continue them for 
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Lemma 4.3. Let b > in the definition of D (4.12 ) be fixed. The minimum ofUo{z) on d^D is achieved at 
n points 

1, j = h 



C — ^25 • • • 5 Cn)' 

The minimal value ofUo{z) on dD is 



-1, j^h 



U \— min Uc\{z) — -. 
zedD ^ 3 



(4.21) 



(4.22) 



Proof. Denote 

d^D := dDf]{z = {zi,Z2,. 
drD := dDf]{z = (zi,Z2^. 

mddiD = drD\JdrD,ie [n]. 

Consider the restriction of Uo{z) on d^D 



: Zi 1}, 



: Zi 



-2-6}, 



n-l 



My) :=t/o((l,y)) = - + 5^F(^.), 



where y = (l/i, 1/2, • • • , Vn-i)- The gradient of t/o is equal to 



ay 



t/o(y) = f (y) := (/(yi), /(2/2), . . . , /(yn-l))^. 



(4.23) 



The definition of / (4.11 ) implies that Uq restricted to d^D has a unique critical point at z = (1, — In-i) 
and Uo{{l^ — In-i)) = 4/3. On the other hand, on the boundary of d^D, dd^D, the minimum of Uo{z) 
satisfies 

4 

min U{z) > -, 
zedd^D 3 



for any 6 > in (4.12). 



Likewise, as follows from the definitions of F (4.16) and D (4.12), for z G d^D, U{z) > 4/3, for 
any choice of 6 > in (4.12). Thus, z = (1, — In-i) minimizes Uq over diD. The lemma is proved by 
repeating the above argument for the remaining faces diD, i e [n] \ {1}. 
□ 



Theorem 4.4. Suppose b > in (4.12) is sufficiently large. There exists 70 > such that for \j\ < Jo, on 
each face d^D^i G [n], U^{z) achieves minimum 

where (p^ : [—70, 70] '^t^^ ^ ^ N' ^ smooth function such that 

d 



Ln-1, 



7=0 



(4.24) 
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and r G ^ is the ith column of the graph Laplacian L after deleting the ith entry. The equations in 
(4.24) are written using the following local coordinates for d^D 



(yi,2/2, . . . ^ (yi,2/2,...,2/z-l,l,yi,...,2/n-2,2/n-l) ^ 9^ D C W. 

Moreover, the minimal value ofU^ on d^D is given by 



:= min t/^ = - + 7 deg{vi) + 0(7^). 



Consequently, 



u, 



'1 



min ?7y = - + 7 min deg{vk) + 0(7^). 



z^dD 



(4.25) 



(4.26) 



Proof Let Uj{y) := Uj{{l, y)), V ^ K"""^ denote the restriction of on d^D: 

n-l 



U,{y) = UHz{y),Hz{y)) + 5^F(y,) + -, 



(4.27) 



i=l 



where ^= (7/1, ?/2, . . . , z{y) := (1, ?/2, . . . , 

Next, we compute the gradient of U^: 



(4.28) 



where f{y) = (/(yi), /(y2), • • • , fiVn-i))- Further, 



^2/ 



;i?z(y),i?z(y)) = 2 



1 
\ 





( ^ \ 




yi 


L 


y2 


J 


\ yn-i J 



2{L^y + l'), (4.29) 



where L^,i e [n], stands for the matrix obtained from L by deleting the ith row and ith column. By plugging 

d 



(4.29»in(4.28), we have 



dy 



U,{y) = j{L'y + l')-f{y). 



The equation for the critical points has the following form 

i?(7,y) :=^{L'y + l')-f{y)=0. 

Note that 

R{0, = 0, — i?(7, y) |^=o,y=-l„_i = -i'ln-i + = 21^ ^ 

Above we used the relation 

— L^ln-i = l^, 



(4.30) 

(4.31) 
(4.32) 



(4.33) 



'Note that P =^ 0,i E [n], as long as Vi is not an isolated node of G. In particular, if G is connected then T ^ Vi e [n]. 
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which follows from the fact that the rows of L sum to zero. 



By the Implicit Function Theorem, for small |7|, the unique solution of (4.31 ) is given by the smooth 
function 6^ : [— 7n, 7nl M^~^ that satisfies 



</.Ho) = -i„_i, A^i(7)|^=o = - 

By taking into account, 



-Q:yR{i^y) |7=o,3y=-i„-i 



9 M 

Q;^R{l,y) |7=0,3y=-l„-i 



d d 

— i?(7, y) |^=o,j/=-i„_i = 2^^ and — i?(7, y) |^=o,2/=-i„_i = 2/„_i, 



from ( |4.34| ) we have 



This shows (4.24). To show (4.25), we use the Taylor expansion of Uy. 



t/o(-l„-i)+7 



(4.34) 



(4.35) 



+ 0(7^ 



J 7=0,y=-ln-i 



^ + -ln-i),H{l, -l„_i)) + 0(7^) = ^ + 7 deg{vi) + 0(7'). (4.36) 



By choosing 6 > in (4.12) large enough one can ensure that Uj{(f)^{j)) for 7 G [—70,70], remains 
smaller than the values of U{z) on the boundary z G d^D. To complete the proof, one only needs to apply 
the same argument to all other faces of d^D, and note that on d~D, the values of U^^ 7 G [—70, 7o], can be 



made arbitrarily large by choosing sufficiently large 6 > in ( |4.12| ) 

□ 



Remark 4.5. The second equation in ( |4.24| ) shows that the minima of the potential function lying on the faces 
corresponding to connected cells move towards the common boundaries of these faces, under the variation 
of 7 > 0. 



4.4 The strong coupling regime 



For small I7I, the minima of Uj{z) are located near the minima of the potential function (cf. (4.16)) 



In this subsection, we show that for larger I7I, the minima of Uj{z) are strongly influenced by the quadratic 
term {Hz, Hz), which corresponds to the coupling operator in the differential equation model (4.10). To 
study the minimization problem for I7I 1, we rewrite Uj{z) as follows: 



Thus, the problem of minimizing for 7 ^ 1 becomes the minimization problem for 

U^{z) {Hz, Hz) + A$(2) ^ min, z G d+D, |A| < 1. 



(4.37) 



(4.38) 
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Lemma 4.6. U^{z) attains the global minimum on d^D at z — 1^: 

:= t/°(ln) = 0. (4.39) 

Proof. U^{z) = {Hz^ Hz) is nonnegative, moreover, 

{Hz, Hz) = if and only if z G ker H = spanjln}. 

Finally, In = ker Hf]d^D. 

□ 

Theorem 4.7. Let Ai(L^), z G [n], denote the smallest eigenvalue of matrix U, obtained from L by deleting 
the ith row and the ith column. Then for < A < min^^[^] Ai(L^), U^{z) achieves its minimal value 
on dD at z — 1^: 

f/\ln) = (4.40) 
provided b > in the definition of D (cf. \4.12\ ) is sufficiently large. 

Remark 4.8. By the interlacing eigenvalues theorem (cf. Theorem 4.3.8, GH), Ai(L^) < A2(L), Vz G [n]. 
(Note that U is not a graph Laplacian.) Furthermore, Ai(L^) > 0, Vz G [n], because G is connected 
(cf. Theorem 6.3, |6|). With these observations. Theorem [477] yields an estimate for the onset of synchrony 
in terms of the eigenvalues of L: 

7 > 2(minAi(L^))-^ > 2(A2(L))-^ (4.41) 



Note that ( 4.41[ ) yields smaller lower bounds for the onset of synchrony for graphs with larger algebraic 



connectivity. In particular, for the families of expanders {Gn} (cf. Example [23]), it provides bounds on the 
coupling strength guaranteeing synchronization that are uniform in n G N. 



For the proof of Theorem 4.7 we need the following auxiliary lemma. 



Lemma 4.9. For z ^ d^D there exists i G [n] such that 



U\z) = ly'Vy + A I ^ + X5 i^(l - [ , 



(4.42) 



i=i 



where y = (yi, y2, ■ • • , Un-iY is defined by 



I J G [n-l]\[z-l]. 

Proof. Since z = (^i, ^2, • • • , ^n)^ ^ d^D, Zi = 1 for some z G [n]. Without loss of generality, let 
zi = 1. Then 

^ = ( . ^ ^' ^ = (^1' 2/2, ... , ^n-i)^, < < 3 + 6, j G [n - 1], 

\ J-n-i y J 
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and 



1 A \ n—1 

U\z) = U\y) -Q{y) + + A F(l - y,) 



where the quadratic function Q{y) = {Hy^ Hy). Differentiating Q{y) yields: 



/ 1 
1 



\ 




Ly^2L^y and ^Q{y)^2L\ 



Therefore, 

□ 



\ ... 1 / 

Q{y) = i/T^iy. 



Proof. (Theorem 4.7 ) Let z G By Lemma 4.9 for some i G [n] and y defined in (4.43 ), we have 

U\z) = \y''Vy + A 1 1 + 1^ F(l - y,) I . (4.44) 
We will use the following observations: 



(a) U is a positive definite matrix, and, therefore, 

y'^Vy > Ai(L^)/y. 



(b) For ^ > 0, 



(c) 



4 
3 



3-3 



^(In) 



A4n 



Using (a) and (b), from ( |4.44| >, we have 

- t/\l„) > {2-^\i{V) - A) > 

provided A < 2-^\x{U). 

This shows that z = !„ minimizes U\ on for A < 2~^ minj^r^i Ai(L*). On the other hand, on 



d D, U\ can be made arbitrarily large for any A > provided 6 > in \A.\2\ is sufficiently large. 

□ 
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4.5 Intermediate coupling strength: formation of clusters 



In this subsection, we develop a geometric interpretation of the spontaneous dynamics of ( |2.6| ) and ( |2.7| ). 
After introducing certain auxiliary notation, we discuss how the spatial location of the minima of Uj(z) on 



the surface of the n— cube encodes the most likely activity patterns of ( |2.6| ) and ( |2.7| ). Then we proceed to 
derive a lower bound on the coupling strength necessary for the development of coherent structures. 



Let G [n], 1 < ii < 12 < • • • < ik ^ ^ and define a (n — A:) -dimensional face of D as 



Qn—k 



The union of all (n — A:) -dimensional faces is denoted by 



u 



{n—k 



D. 



I<ii<i2<-"<i/c<n 



(4.45) 



(4.46) 



Suppose ^ is a point of minimum of Uj{z) on d^D. If ^ G 5^ ^D, k > 2, then with high probability 



the network discharges in A:— clusters. The analysis in [43 shows that for small I7I there is practically no 
correlation between the activity of distinct neurons. On the other hand, when the coupling is strong, all 



cells fire in unison (cf. §4.4[ ). Lemma [4^101 provides a lower bound on the coupling strength needed for the 
formation of clusters. 

Lemma 4.10. Let ^ G dD be a point of global minimum ofU^ on dD. If^ G d^^~^^D for some k>2 then 

2 



7> 



(4.47) 



Proof Suppose z e ^D. Without loss of generality, we assume that 

z^ln- (0, 0, ^2, 2/3, ... , Vn-i) In " (0, 0, y), ^ G M"-' > 0, z G [n - 1] \ {1}. 
Denote y = (yi, ^2, 1/3, ... , Vn-i) = (l/i, y), Vi > 0, and z = In-i - (0, y). Thus, 

n—l . 



U^{z) = ^y'"Liy + ^F(l-yi) + 

i=l 



U,{z) 



n-1 



1=2 



(4.48) 
(4.49) 



where L^^ is a matrix obtained from L by deleting the first and the second rows and columns. The Laplacian 
of G can be represented as 

L = diag(deg(i;i),deg(i;2), . . . ,deg(i;^)) - A, 



where the adjacency matrix A is nonnegative (cf. (2.1 1 )). Therefore, for nonnegative y ^ {yi^y) ^ 



y L y-y L y < deg {v2)y2 < max deg {vk)y2, 

keln] 



(4.50) 
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Further, for any < S < 1, 

n-l 



1=1 K i=2 ) 

provided < y2 < 3(5. The combination of ( |4.48| )-( |43T] ) yields 



fory2 e (0,35). By (|432), 
provided 



min U^iz) > min U^iz), 



7 < 



2(1 -<5) 



(4.51) 



U^iz) < U^{z) + J max deg (vk) -yl + \yl < U^{z) + ( J max deg {vk) - 1 + s] yl (4.52) 

Z ke[n\ o \Z ke[n\ J 



(4.53) 



max;,^[^] deg [v].)' 

The statement of the lemma follows from the observation above by noting that 5 G (0, 1) in (4.53) is 
arbitrary. 

□ 



5 An alternative view on synchrony 



The variational analysis of the previous section shows that when the coupling is strong (7 > (A2(L)) the 



neurons in the network fire together (cf. Theorem 4.7 ). In this section, we use a complementary approach to 
studying synchrony in the coupled system. We show that strong coupling brings about the separation of the 
timescales in the system's dynamics, which defines two principal modes in the strong coupling regime: the 
fast synchronization and slow large-deviation type escape from the potential well. The analysis of the fast 
subsystem elucidates the stability of the synchronization subspace. In particular, it reveals the contribution 
of the network topology to the robustness of synchrony. The analysis of the slow subsystem yields the 
asymptotic rate of the network activity in the strong coupling regime (cf. Theorem [4?7|). 



5.1 The slow-fast decomposition 



Our first goal is to show that when coupling is strong the dynamics of the coupled system ( |4.10| ) has two 
disparate timescales. To this end, we introduce the following coordinate transformation 



pn— 1 



where 



and U G 



^ Rz and 77 = n ^In^^, 



(5.1) 



(n— l)xn 



is the coboundary matrix corresponding to the spanning tree G of G (see (2.9)) 
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Lemma 5.1. Equation ( [5.7p defines an invertible linear transformation: 

Z = 7/ln + Si. 

Matrix S = {sij) G R^><(^~1) satisfies 

\sij\ < 1, G [nf and IxJ S — 0. 



(5.2) 



(5.3) 



Proof. Fix z G [n]. For each j G [n] \ {z} there exists a unique path connecting nodes vi G V^(G) and 



Vj G y (G) and belonging to G 



n-l 



p{hj) = ^(^k{hjyk, ^ {o,±i}. 



Thus, 



k=i 



n-l 



(5.4) 



(5.5) 



k=i 



where ^ = (Ci, ^2, • • • , Cn-i)- By summing n — 1 equations (5.5 ), adding the identity Zi = Zi to the resultant 
equation, and dividing the result by n, we obtain 



n—l n—1 

= ^ + X^^i/cC/c, where Sik = -n'^'^akiij) 

k=l 3=1 



(5.6) 



The first inequality in (5.3) follows from the formula for in(5.6)and |cr/e(i, j)| < 1. To show the second 
identity in (|5.3|), add up equations (|5.6|) for i G [n] and use the definition of 77: 



□ 



Lemma 5.2. Suppose G is a connected graph and L G W^'^^ is its Laplacian. There exists a unique 
L G r(^-i)x(^-i) such that 

HL = LH, (5.7) 

where H G r(^~^)x^ is the coboundary matrix of the G d G, a spanning tree of G. The spectrum of L 
consists of all nonzero eigenvalues of L 



spec(L) = spec(L) \ {0}. 



(5.8) 



Proof. Since G is connected and G is a spanning tree of G, 

rank H — n — 1 and ker H — ker L — span {In}- 



(5.9) 



The existence and uniqueness of the solution of the matrix equation (5.7 ), L, is shown in Lemma 2.3 of ll29l . 



Equation follows from Lemma 2.5 of ||29j 

□ 



We are now in a position to rewrite (4.10) in terms of r;). 
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Lemma 5.3. The following system of stochastic dijferential equations is equivalent in law to (4.10): 



7 



n 



where 

Qi(e,?/) = 2r;e + diag(6,6,---,en-i)5C, and Q2{i) ^ rT^C Si, 

/ Rowi(5) + Row2(5) \ 
Row2(5) + Row3(5) 

VRow„_i(5) + Row,(5)/ 
Throughout this section, W and w denote the white noise processes in R" and R respectively. 



(5.10) 
(5.11) 

(5.12) 



Proof. After multiplying both sides of (4.10> by H and using (5.7 1, we have 



i = -^-U + m{z) + aHWt, 



(5.13) 



where 



m{z) 



( 4-4 \ 

\4i ~ 4i-l/ 



I {zi + Z2)ii \ I (2r?+[Rowi(5) + Row2(5)]0a \ 
{Z2 + ^3)6 ^ (2r? + [Row2(5) + Row3(5)] 06 

\{Zn-l + Zn)in-J \(2?? + [Row„_i(5) + Row„(5)] O^n-lJ 



This shows ( 5.10 >. By multiplying ( 4.10 > by n ^l^J, using !„ G ker L^, and the definition of rj, we have 

n 

i=l 

Here, we are using the fact that the distributions of n~^lyJW and ~ -^w coincide. Using the definition 
of / ^J\\ and (|53]), we have 



J2 {fiv + Row,(^)0 - fiv)} = 2r?(l„"^5e) + Tt (SOiSO'' = fS^S^. 
i=i 



□ 



5.2 Fast dynamics: synchronization 



For 7 » 1, the system of equations ( |5.10| ) and ( |5.11| > has two disparate timescales. The stable fixed point 
^ = of the fast subsystem (|5.10[> corresponds to the synchronous state of (|4.1Q[> 



Zl = Z2 
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In this section, we analyze the stability of the steady state of the fast subsystem. Specifically, we identify 
the network properties that determine the rate of convergence of the trajectories of the fast subsystem to the 
steady state and its degree of stability to random perturbations. These results elucidate the contribution of 
the network topology to the synchronization properties of the coupled system ( |4.10| ). 

By switching to the fast time ( |5.10| ) and ( |5.11| ), we have 

X = -LX + 5gi(x,y) + aW, 

Y 



\/n 



(5.14) 
(5.15) 



where 



X{s) = C{Ss), Y{s) = ri{Ss), < 5 = 27"^ < 1. (5.16) 
The leading order approximation of the fast equation ( |5.14| ) does not depend on the slow variable Y: 



X 



-LX + dHW. 



(5.17) 



The solution of (5.17 ) with deterministic initial condition X(0) = x G ^ is a Gaussian random process. 

(5.18) 



The mean vector and the covariance matrix functions 

m{s) := E X{s) and V{s) := E \{X{s) - m{s)){X{s) - m{s)y 
satisfy linear equations 



m^-Lm and V ^ LV ^VL ^- d^HH'^ (5.19) 

Recall that L is a positive definite matrix, whose smallest eigenvalue Ai(L) is equal to the algebraic con- 
nectivity of G, a = \2{L) (see Lemma 5.2). By integrating the first equation in (5.19), we find that 

E X[s) = exp{-5L}x < Ci expj-asjx ^0, 5 ^ oc, (5.20) 

for some Ci > 0. Thus, the trajectories of the fast subsystem conv erge in mean to the stable fixed point 
X = 0, which corresponds to the synchronization subspace of (4.10). Moreover, the rate of convergence is 
set by the algebraic connectivity a. 



Further, 



n-l 



TtV{s) = E\X{s) -m{s) 



^var X(5) + o(l), 5 > 1, 

i=l 



(5.21) 



measures the spread of the trajectories around the synchronization subspace. By integrating the second 
equation in ( |5.19| ), we have 



a-^TTV{s) = Tr 
= Tr 



/ exp{{u — s)L}Aexp{{u — s)L}d' 
Jo 

A / exp{—2L}udu 
. Jo 



where 



/^(G, G) := Tr {^-^A} and A = HH^. 



(5.22) 



(5.23) 
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Parameter G) quantifies the mean square stability of the synchronization subspace. In ^5.4 
that G) depends on the properties of the cycle subspace of G. 



we show 



For small a > 0, a typical trajectory of ( |5.17| ) converges to a small neighborhood of the stable equilib- 
rium at the origin. However, eventually it leaves the neighborhood of the origin under persistent random 
perturbations. Next, we estimate the time that the trajectory of the fast subsystem spends near the origin. 



Let p > and Bp = {X e 



pn—l 



\X\ < p], p > 0. By X{t) we denote the solution of (5.17) 



satisfying initial condition X(0) = x G Bp. Define the first exit time of the trajectory of (5.17) from Bp\ 

T(X,p) = inf{X(t) = p}. (5.24) 

Using the large deviation estimates (cf. HH), we have 



where 



lim P ^ {exp{a-2(T/o(p) - h)} < t(x, p) < exp{a-\Vo(p) + h)}} = 1, V/i > 0, 



Vo{p) = min ^{Lx,x) = ^ 

\x\=p Z Z 



(5.25) 



(5.26) 



is the minimum of the positive definite quadratic form 2 ^ (Lx, x) on the boundary of Bp. 



By combining (5.25 ) and (5.26) we obtain the logarithmic asymptotics of the first exit time from Bp 

t(X,p) X exp 



ap 



a2 



(5.27) 



For small S > (i.e., for large 7 » 1), ( |5.27| ) applies to the fast equation ( |5.14| ). Switching back to the 

(5.28) 



original time, we rewrite the estimate for the first exit time for the fast equation ( |5.10| ): 

t(^,p) X exp 



-3-\-L 



Finally, choosing p := 7 2 for an arbitrary fixed < ^ < |, we have 



t(C, p) X exp{0(a "Y)}, 7 > 1, < a < 1. 



(5.29) 



5.3 The slow dynamics: escape from the potential well 



We recap the results of the analysis of the fast subsystem. The trajectories of the fast subsystem ( |5.10| ) enter 

— 3+t- 

an 0(7^~) neighborhood of the stable equilibrium ^ = (corresponding to the synchronization subspace 
of the full system) in time 0(7"^ In 7) and remain there with overwhelming probability over time intervals 
0{exp{0{a~'^Y)}' During this time, the dynamics is driven by the slow equation (5.11 ), which we analyze 
next. 



While the traj ectory of the fast subsystem stays in the 0(7^"^) neighborhood of the equilibrium, the 
quadratic term in (5.11) Q2{0 ^ O(7~^^0 is small. Thus, the leading order approximation of the slow 
equation is independent of ^ on time intervals 0(exp{0((j~^)}) 



a 



n 



(5.30) 
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Equation (5.30) has a stable fixed point at the origin, where the potential function ^(77) attains its min- 



imum value (see (4.16)). The escape of the trajectories of (5.30) from the potential well, defined by ^(77) 
corresponds to spontaneous synchronized discharge of the coupled subsystem. Suppose 77(0) = 770 < 1 and 
define 

r{r]) = inf{t > : r]{t) = 1}. (5.31) 
By the large deviation theory, we have 



r{r]) 



r 2AFn ] 
exp <^ 



AF:=F(1)-F(0) 



4 
3' 



(5.32) 



Equation ( |5.32| ) provides the estimate for the frequency of spontaneous activity in the strong coupling 
regime. Note that ( |5.32| ) is consistent with the estimate in Theorem |4.7| derived using the variational ar- 
gument. Therefore, the analysis in this section yields the dynamical interpretation for the minimization 



problem for (4.20) in the strong coupling regime. 



We summarize the results of the slow-fast analysis. In the strong coupling regime, the dynamics splits 
into two modes: fast synchronization and slow synchronized fluctuations leading to large deviation type 
discharges of the entire network. By analyzing the fast subsystem, we obtained estimates of stability of the 
synchronization subspace. The analysis of the slow subsystem yields the asymptotic estimate of the firing 
rate in the strong coupling regime. 



5.4 The network topology and synchronization 



The stability analysis in §5.2| provides interesting insights into what structural properties of the network are 
important for synchronization. In this subsection, we discuss the implications of the stability analysis in 
more detail. 



First, we rewrite ( |5.20| ) and ( |5.22| ), using the original time 

E i{t) < Ci exp{-2"^a7"4}x, 



E 



(5.33) 
(5.34) 



where ^{t) = X(5 ^t). Equation (5.33 ) shows that the rate of convergence to the synchronous state depends 



on the coupling strength, 7, and the algebraic connectivity of the network, a. The convergence is faster for 
stronger coupling and larger a. This, in particular, implies that the rate of convergence to synchrony in 



networks on spectral expanders (cf. Example 2.5) remains 0(1) as the size of the network grows without 
bound. In particular, the families of the random graphs (cf . Example |2.4[ ) have this property. In contrast, 
for many networks with symmetries (cf. Example [Z2|) the algebraic connectivity a = o(l)as7T.^oc, 



and, therefore, by ( |5.33| ), synchronization requires longer time, if the size of the network grows. For a more 
detailed discussion of the synchronization properties of the networks on expanders, we refer the interested 
reader to |[29l . 



Next, we turn to Equation ( |5.34| ), which characterizes the dispersion of the trajectories around the syn- 
chronization subspace. E may be viewed as a measure of robustness of synchrony to noise or, more 
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generally, to constantly acting perturbations. By ( |5.34[ ), the synchrony is more robust for stronger coupling, 

2 

because the asymptotic value of E ^{t) ^ as 7 ^ oc. The contribution of the network topology to 

the mean-square stability of the synchronous state is reflected in /^(G, G). Trajectories of the networks with 
smaller G) are more tightly localized around the synchronization subspace. 

To explain the graph-theoretic interpretation of /^(G, G), we review the structure of t he cy cle subspace 
of G. Recall that the edge set of the spanning tree G consists of the first n — 1 edges (see (2.9)): 

E{G) = {ei,e2, . . • ,e^_i}. 

If G is not a tree, then 

E{G) \ E{G) = {en, e^+i, . . . , e^+c-i}, 

where c is the corank of G. To each edge e^+j^ , A: G [c], there corresponds a unique cycle Ok of length \Ok\, 
such that it consists of and the edges from E{G). The following lemma, relates the value of G) to 
the properties of the cycles {Ok}k=i' 

Lemma 5.4. I{29\l Let G = {V{G),E{G)), \V{G)\ = n, be a connected graph. 
A If G is a tree then 

k{G,G) = (5.35) 
B Otherwise, let G C G be a spanning tree ofGa and {Ok}k=i be the corresponding independent cycles. 
B.l Denote 



k=i 



Then 



B.l IfO <c<n-l then 



1 <!^<l. (536) 



where 



n-1 \ M J ~ n-1 
M = max{|Ofc| + VlO^nOH}. 
B.3 IfOk^k G [c] are disjoint. Then 



1 + /X n — 1 



1 - < ^ < 1. (5.37) 



'"°'°'^i-^(i-iEio^|-')- 



n — 1 n 



In particular, 
and 

i^{G,G) 



n{G,G) c 
n-1 - n-lV^ min^^fc] \Ok 



>1 r l-^ ^ >1 



n-1 n - 1 I E/cG[c] |0/c| / n - 1 V max;.^[c] |Oa; 
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The asymptotic estimate of the mean s quar e stabihty of the synchronization subspace in ( |5.34| ) combined 



with the estimates of a^(G, G) in Lemma |5.4| show how the structure of the cycle subspace of G translates 
into the stability of the synchronous state. In |29|, one can also find an estimate of the asymptotic stability 
of the synchronization subspace in terms of the effective resistance of the graph G. These results show how 
the structural properties of the network shape synchronization properties of the coupled system. 



6 Discussion 



In this paper, we presented a variational method, which reduces the problem of pattern-formation in electri- 
cally coupled networks of neurons to the minimization problem for the potential function U^{z) on the sur- 
face of a unit n— cube, dD. The variational problem provides geometric interpretation for the spontaneous 
dynamics generated by the network. Specifically, the location of the points of minima of the constrained 
potential function U^{z) = Uj{z e dD) corresponds to the most likely patterns of network activity. 

The variational formulation has several important implications for the analysis of the dynamical prob- 
lem. First, the minimization problem has an intrinsic bifurcation structure. By a bifurcation of the problem 



( |4.20| ), it is natural to call the value of the parameter 7 = 7* corresponding to the structural changes in the 
singularity set of U^{z). For example, the number of the minima of U^{z) can change due to collisions 
of the singularities with each other or with the boundaries of the faces of a given co-dimension (a border 
collision bifurcation). In either case, the qualitative change in the configuration of the po ints of minima of 



Uj{z) signals the transformation of the attractor of the randomly perturbed system (4.10). The study of the 



constrained minimization problem ( |4.20| ) identified three main regimes in the network dynamics for weak, 
intermediate, and strong coupling. These results hold for any connected network. We expect that under 
certain assumptions on the network topology (e.g., in the presence of symmetries, or alternatively in net- 
works with random connections), a more detailed description of the bifurcation events preceding complete 
synchronization should be possible. 

The analysis of the variational problem has also helped to obtain quantitative estimates for the network 
dynamics. For both weak and strong coupling, we derived the asymptotic formulae for the dependence of 
the firing rate as a function of the coupling strength (see ( |4.26| ) and ( |4.40[ )). Surprisingly, in each of these 



cases, the firing rate does not depend on the structure of the graph of the network beyond its degree and 
order. In particular, the network of equal degree with connectivity patterns as different as symmetric and 
random exhibit the same activity rate (see Fig.lSl 



The geometric interpretation of the spontaneous dynamics yields a novel mechanism of the formation of 
clusters. It shows that the network fires in A:— clusters, whenever Uj{z) has a minimum on a co-dimension 
k G [n] face of dD. In particular, the network becomes completely synchronized, when the minimum of 
Uj{z) reaches In G d^D n ker(i7). This observation allows one to estimate the onset of synchronization 



(cf. Theorem 4.7 ) and cluster formation (cf. Lemma 4. 10). Furthermore, we show that in the strong coupling 



regime, the network dynamics has two disparate timescales: fast synchronization is followed by an ultra- 
slow escape from the potential well. The analysis of the slow-fast system yields estimates of stability of the 
synchronous state in terms of the coupling strength and structural properties of the network. In particular, it 
shows the contribution of the network topology to the synchronization properties of the network. 

We end this paper with a few concluding remarks about the implications of this work for the LC net- 
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Figure 1 1 : The responses to stimulation of the two networks at different values of the coupling strength: 
^ = (left column) and g = 0.02 (right column). The networks shown in a and b generate different spatio- 
temporal patterns. However, the firing rates corresponding to these activity patterns are close (see c and d). 
When 10 neurons in the middle of each of these networks receive a current pulse, the network responses 
differ (see e and f). The firing rate in the first network during the stimulation changes very little (see g), 
while the response of the second network is clearly seen from the firing rate plot (see g). 
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work. The analysis of the conductance-based model of the LC network in this paper agrees with the study of 
the integrate-and-fire neuron network in ll38]| and confirms that the assumptions of spontaneously active LC 
neurons coupled electrically with a variable coupling strength are consistent with the experimental obser- 
vations of the LC network. Following the observations in |38 | that stronger coupling slows down network 
activity, we have studied how the firing rate depends on the coupling strength. We show that strong coupling 
results in synchronization and significantly decreases the firing rate (see also |[32l [34l)- Surprisingly, we 
found that the rate can be effectively controlled by the strength of interactions already for very weak cou- 
pling. We show that the dependence of the firing rate on the strength of coupling is nonmonotone. This has 
an important implication for the interpretation of the experimental data. Because two distinct firing patterns 
can have similar firing rates, the firing rate alone does not determine the response of the network to exter- 
nal stimulation. This situation is illustrated in Fig. [TT] We choose parameters such that two networks, the 
spontaneously active (Fig. [TT^ ) and the nearly synchronous one (Fig. [TTJ3), exhibit about the same activity 
rates (see Fig. [TT]:,d). However, because the activity patterns generated by these networks are different, 
so are their responses to stimulation (Fig. [TT^,f). The network in the spontaneous firing regime produces 



a barely noticeable response (Fig. [TTg), whereas the response of the synchronized network is pronounced 
(Fig. (TTJi). Network responses similar to these were observed experimentally and are associated with the 
good (Fig. [TT]i) and poor (Fig. [TT^) cognitive performance |[38l . Our analysis suggests that the state of the 
network (i.e., the spatio-temporal dynamics) rather than the firing rate, determines the response of the LC 
network to afferent stimulation. 



The main hypotheses used in our analysis are that the local dynamical systems satisfy Assumption [4T 
and interact through electrical coupling. The latter means that the coupling is realized through one of the 
local variables, interpreted as voltage, and is subject to the two Kirchhoff 's laws for electrical circuits. In this 
form our assumptions cover many biological, physical, and technological problems, including power grids, 
sensor and communication networks, and consensus protocols for coordination of autonomous agents (see 
|[29ll and references therein). Therefore, the results of this work elucidate the principles of pattern formation 
in an important class of problems. 

Acknowledgments. This work was partly supported by the NSF Award DMS 1 109367 (to GM). 



Appendix. The parameter values used in the biophysical model ( |2J] ) and ( |2^ 



To emphasize that the results of this study do not rely on any specific features of the LC neuron model, 
in our numerical experiments we used the Morris-Lecar model, a common Type I biophysical model of an 
excitable cell ll35i . This model is based on the Hodgkin-Huxley paradigm. The function on the right hand 
side of the voltage equation ( |2.1| ), lion = Ica + Ik + Ih models the combined effect of the calcium and 
sodium currents, Ica, the potassium current, Ik, and a small leak current, 

ICa{v) gCamoo{v){v - Eca)^ 

Ixiv^n) = gKn{v-EK), 
Ii{v) = gi{v-Ei). 



Constants Eca, Ek, and Ei stand for the reversal potentials and gca, 9k, and gi denote the maximal 
conductances of the corresponding ionic currents. The activation of the calcium and potassium channels are 



32 



modeled using the steady-state functions 

— Ul 



^oo('^) = 0-5 1 + tanh 



^2 



and noo{v) = 0.5 ( 1 + tanh 



V - 1^3 
U4 



and the voltage-dependent time constant 

t{v) = ^cosh 

The parameter values are summarized the following table. 

Table 



Eca 


llOmV 


9K 


Ss-' 


Er 


-UmV 


1^2 




gca 


4s-i 


</> 


0.067 


m 


2s-i 


El 


-60mV 


c 


20 iiF/crn^ 




-l.2mV 




l2mV 


U4 


llAmV 



33 



References 

[1] G. Aston- Jones and J.D. Cohen, An integrative theory of Locus Coeruleus-Norepinephrine function: 
adaptive gain and optimal performance, Annu. Rev. NeuroscL, 28:403-50, 2005. 

[2] V. Alvarez, C. Chow, E.J. Van Bockstaele, and J.T. Williams, Frequency-dependent synchrony in locus 
coeruleus: role of electronic coupling, Proc. Nat. Acad. Sci. USA, 99, 4032-4036, 2002. 

[3] N. Berglund and B. Gentz, Noise-Induced Phenomena in Slow-Fast Dynamical Systems: A Sample- 
Paths Approach, Springer, 2006. 

[4] Berridge, C.W. and Waterhouse, B.D., The locus coeruleus-noradrenergic system: modulation of be- 
havioral state and state-dependent cognitive processes. Brain Research Reviews, 42, 33-84, 2003. 

[5] Brown, Eric and Moehlis, Jeff and Holmes, Philip and Clayton, Ed and Rajkowski, Janusz and Aston- 
Jones, Gary, The influence of spike rate and stimulus duration on noradrenergic neurons. Journal of 
Comp. Neuroscience,\l , 13-29, 2004. 

[6] N. Biggs, Algebraic Graph Theory, second ed., Cambridge University Press, 1993. 

[7] Bela BoUobas, Modern graph theory. Graduate Texts in Mathematics, vol. 184, Springer, New York, 
1998. 

[8] S.-N. Chow and J.K. Hale, Methods Of Bifurcation Theory, Springer- Verlag New York Inc, New York, 
1982. 

[9] F.R.K. Chung, Spectral Graph Theory, CBMS Regional Conference Series in Mathematics, No. 92, 
1997. 

[10] B.W. Connors and M.A. Long, Electrical synapses in the mammalian brain, Annu. Rev. NeuroscL, 
27:393-418, 2004. 

[11] S. Coombes, Neuronal networks with gap junctions: A study of piece-wise linear planar neuron mod- 
els, SIAM Journal on Applied Dynamical Systems, 7, 1101-1129, 2008. 

[12] Martin V. Day, On the exponential exit law in the small parameter exit problem, Stochastic s^ 8, 297- 
323, 1983. 

[13] M.I. Freidlin, On stable oscillations and equilibriums induced by small noise, /. of Stat. Phys., 103 
(1-2), 283-300, 2001. 

[14] M.I. Freidlin and A.D. Wentzell, Random perturbations of dynamical systems, 2nd ed., Springer, New 
York, 1998. 

[15] J. Friedman, A Proof of Alon's Second Eigenvalue Conjecture and Related Problems, Memoirs of the 
American Mathematical Society, vol. 195, 2008. 

[16] M. Fiedler, Algebraic connectivity of graphs. Czech. Math. J. 23(98), 1973. 

[17] J. Jost, Dynamical networks, in J. Feng, J. Jost, and M. Qian, (Eds.) Networks: From Biology to Theory, 
Springer, 2007 



34 



[18] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vec- 
tor Fields, Springer, 1983. 

[19] Juan Gao and Philip Holmes, On the dynamics of electrically-coupled neurons with inhibitory 
synapses, /. Comput. NeuroscL, 22:3961, 2007. 

[20] J.K Hale, Ordinary Differential Equations, Krieger, 2nd edition, 1980. 

[21] R.A. Horn and C.R. Johnson, Matrix Analysis, Cambridge University Press, 1999. 

[22] S. Hoory, N. Linial, and A. Wigderson, Expander graphs and their applications. Bulletin of the Ameri- 
can Mathematical Society, vol. 43, no. 4, pp 439561, 2006. 

[23] I. Karatzas and S.E. Shreve, Brownian Motion and Stochastic Calculus, 2nd ed.. Springer, New York, 
1991. 

[24] James P. Keener, Propagation and its failure in coupled systems of discrete excitable cells, SIAM Jour- 
nal on Applied Mathematics, v.47 n.3, p.556-572, 1987. 

[25] Yuri A. Kuznetsov, Elements of applied bifurcation theory. Springer, 1998. 

[26] T. Lewis and J. Rinzel, Dynamics of spiking neurons connected by both inhibitory and electrical cou- 
pling, /. Comp. NeuroscL, 14:283-309, 2003. 

[27] A. Lubotzky, R. Phillips, and P. Sarnak, Ramanujan graphs, Combinatorica, 8, 161-278, 1988. 

[28] G. Margulis, Explicit group-theoretic constructions of combinatorial schemes and their applications in 
the construction of expanders and concentrators. (Russian) Problemy Peredachi Informatsii 24 (1988), 
no. 1, 51-60; (English translation in Problems Inform. Transmission 24 (1988), no. 1, 39^6). 

[29] G.S. Medvedev, Stochastic stability of continuous time consensus protocols, submitted, arXiv preprint: 
1007.1234. 

[30] G.S. Medvedev, Synchronization of coupled limit cycles, /. Nonlin. ScL, 21, 441^64, 201 1. 

[31] G.S. Medvedev, Synchronization of coupled stochastic limit cycle oscillators. Physics Letters A (374), 
1712-1720, 2010. 

[32] G.S. Medvedev, Electrical coupling promotes fidelity of responses in the networks of model neurons. 
Neural Computation, 21 (11), 3057-3078, 2009. 

[33] G.S. Medvedev and N. Kopell, Synchronization and transient dynamics in the chains of electrically 
coupled FitzHugh-Nagumo oscillators, SIAM J. ofAppl Math., 61, 1762-1801, 2001. 

[34] G.S. Medvedev and S. Zhuravytska, Shaping bursting by electrical coupling and noise, submitted, 
arXiv preprint: 1111 .0642. 

[35] J. Rinzel and G.B. Ermentrout, Analysis of neural excitability and oscillations, in C. Koch and I. Segev, 
eds Methods in Neuronal Modeling, MIT Press, Cambridge, MA, 1989. 

[36] P. Sarnak, What is an expander?. Notices of the American Mathematical Society, 51, 762-763, 2004. 



35 



[37] S.J. Sara, The locus coeruleus and noradrenergic modulation of cognition, Nat. Rev. NeuroscL, 10, 
211-223, 2009. 

[38] M. Usher, J.D. Cohen, D. Servan-Schreiber, J. Rajkowski and G. Aston- Jones, The role of locus 
coeruleus in the regulation of cognitive performance. Science, 283 (1999), pp. 549-554. 



36 



